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Abstract 

A relaxation method is employed to study a rotating dense Bose-Einstein conden- 
sate beyond Thomas-Fermi approximation. We use a slave-boson model to describe 
the strongly interacting condensate and derive a generalized non-linear Schrodinger 
equation with kinetic term for the rotating condensate. In comparison with previ- 
ous calculations, based on Thomas-Fermi approximation, significant improvements 
are found in regions, where the condensate in a trap potential is not smooth. The 
critical angular velocity of the vortex formation is higher than in the Thomas-Fermi 
prediction. 

1 Introduction 

The appearance of vortices in a rotating dilute Bose-Einstein condensate (BEC) is a phe- 
nomenon that is intrinsically related to the existence of superfluidity [HE]. This has been 
intensively studied j3j H] after the achievement of BEC in alkali-metal gases [5j |6] . The 
vortex formation can be understood as a local suppression of the BEC due to rotation. It 
happens when the angular velocity exceeds a critical value Q c [7j. The critical angular 
velocity is obtained by an energy argument [7J, in which the energy of the BEC without 
a vortex is higher than the energy of the BEC with a vortex for O > £l c . With higher 
angular velocities more than one vortex can be created to minimize the energy of the BEC. 
In this paper we will concentrate on values of O where only the the creation of a single 
vortex is possible. 

Since vortices are accompanied by a local suppression of the BEC, other mechanisms 
of BEC suppression may interfere with the vortex formation. An important example is 
the depletion of the BEC by a strong interaction between the bosonic particles. The latter 
is possible in the form of repulsive s-wave scattering in a dense BEC 

In the case of a dilute BEC (i.e. for weak interaction) a mean-field theory in terms 
of the Gross-Pitaevskii equation provides a description of the vortex formation [8]. The 
relation between the number of particles in the BEC and the critical angular velocity was 
determined within this approach [9j. This can be extended to a dense Bose gas, where not 
all particles contribute to the BEC (i.e. there is a depletion of the BEC). The discussion 
of the dense Bose gas requires a statistical description of the entire Bose gas, using, 
for instance, a functional-integral representation [10]. Nevertheless, it turns out that a 
classical field equation for the BEC order parameter can also be derived, which is equivalent 
with the Gross-Pitaevskii equation in the limit of a dilute Bose gas [11] and a renormalized 
Gross-Pitaevskii equation for higher densities. In the regime of high densities the strongly 
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interacting Bose gas had been studied within this statistical approach [TOJ, [III EJ [TOJ [TU 
[TOJ,[TO]. A relatively simple method to study strong interaction is the slave-boson approach. 
This method was originally suggested for strongly interacting fermions [17\ [18j and was 
later also applied to bosons [10]. The central idea of the method is to construct operators 
or fields for particles and holes in a fictitious vacuum, and to introduce a constraint, 
which guarantees that each position is occupied either by a particle or by a hole. This 
leads to a dynamics, where tunneling of physical particles is described by an exchange 
of a particle with a hole, and no multiple occupation by particles is possible. Although 
this construction makes the theory more complex (because of the additional holes), it 
also provides more freedom to choose a proper mean-field approximation. Moreover, the 
bosonic version of the slave-boson approach, which will be applied in this paper, allows 
for an exact treatment of the constraint [TO]. Here we will not present the details of the 
slave-boson approach but refer to the literature. A recent discussion can be found in Ref. 
|15| . Within a mean-field (Thomas-Fermi) approximation of a slave-boson theory we were 
able to study the formation of a vortex in the BEC of a dense Bose gas [15] , However, a 
mean-field approach is a reliable approximation only for a BEC that is smoothly changing 
in space. The latter is not the case near the core of a vortex, where the condensate density 
is changing strongly on short scales. Therefore, it is important to extend the mean-field 
theory to a classical field theory with a kinetic term. This leads to a nonlinear partial 
differential equation which can be treated with a numerical method. In the present work 
we are using a relaxation approach. Relaxation techniques have been extensively used 
to obtain the ground-state solution of the Gross-Pitaevskii equation [211 EH E31 [Ml EH] • 
These techniques require interaction and need little computational efforts. The method 
provides stable solutions of the slave-boson approach. 

This work is organized as follows. In Section 2 we summarize the main equation 
of slave-boson approach. The trapped condensate is shown in spherical and cylindrical 
coordinates in Section 3. In Section 4 we discuss the numerical procedure used in this 
work and present the numerical results with and without vortex. Afterwards we review 
our main results in section 5. 



2 Slave-Boson Approach 



We consider a grand-canonical ensemble of bosons with chemical potential fi at tempera- 
ture T = l/ksP- Its thermodynamics is described by the partition function 
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where H is the Hamiltonian of the interacting Bose gas, N the particle number operator, 
and Tr is the trace over all possible quantum states of the system. If the interaction is due 
to a hard core of a given radius a, the Bose gas can be approximated by a lattice gas with 
lattice constant a. Then we can apply the slave-boson approach [TOJ CEH El EE HH [TOJ [TO] 
which allows us to write the partition function Z as a functional integral with respect to 
a complex (order-parameter) field $(r). Here we start with an expression for Z, whose 
derivation was given in Ref. [TO] : 
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a is a free parameter which appears as an ambiguity in the functional integral [19]. It is 
necessary for the integration, but the exact functional integral Z does not depend on the 
choice of a. However, approximations may depend on it. In agreement with a previous 
work [11] we choose a = 1/5.5 for the following calculation. In S we have introduced the 
parameters 

= aJp , // = A= , (5) 



where J is the tunneling rate of the lattice model (i. e. before taking the continuum 
limit). 

Since $ is the order-parameter field, the number of condensed bosons can be written 
as an integral over |<3?(r)| 2 : 



Nn 



[ |$(V)| 2 d 3 r 

a* (l + l/a) 2 i 1 1 >l a 



(6) 



Thus we can associate a 3 j^p with the density of condensed bosons. 

The evaluation of the partition function in Eq. ([2D requires a functional integration. 
This can be performed approximately by a saddle-point integration. The saddle point is 
determined as the solution of a vanishing variation of S (i.e. 5S = 0). This leads to 
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where the last term can be explicitly evaluated as 
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Eq. d?]) corresponds to the Gross-Pitaevskii equation 
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The latter is obtained as a special case of Eq. ([7|) for a weakly interacting Bose gas 
|15j . Here M is the mass of the bosons, and the interaction strength of the two-particle 
interaction g is related to the s-wave scattering length a s by 
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There is a direct correspondence between the parameters J, fi, and a in Eq. ((Jj), and the 
parameters M, jl, and g in Eq. (fTUj) [16J 
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Since Eq. ([7]) describes the condensate for arbitrary density, we will use this in the 
subsequent calculation to evaluate the condensate profile of a trapped and rotating dense 
Bose gas. 
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3 Dense condensate in a spherical harmonic trap potential 



We assume a spherical trap potential given by 



2 2 
^ho r 
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where Who is the trap angular frequency measured in rad/s. In typical experiments, the 
oscillator length d^ = y/h/Moj^o is of the order of a few fim |26j. The external potential 
can now be included in a space-dependent chemical potential as fi — ► fj,(r) = /j, — V(r). 

Considering, for instance, 85 Rb atoms near a Feshbach resonance [27], we can study a 
Bose gas in a dense regime with a scattering length a s ~ a ~ 200nm. In our calculations 
we choose the parameters 
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in order to compare with previous results in Ref. |15j . 

In a spherical symmetric trap we can assume that the condensate is also spherically 
symmetric: 

$(r) = <j>{r) . (15) 

Then Eq. ([7]) takes the form 
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Now we consider a Bose gas which rotates with angular velocity f2 about the z-axis. In 
the rotating frame this can be described by adding —h£lL z to the kinetic energy in Eq. 
(J7J), where the z-component of the angular momentum operator is 
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and <f is the polar angle in cylindrical coordinates. 

For a condensate with a single vortex along the z-axis, we use cylindrical coordinates 
(rj_, z, ip) and make the ansatz 
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for the condensate order parameter. For m = there is no vortex and for m = 1 there is 
a single vortex. Then the classical field equation reads 
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where Q! is the rescaled angular velocity 
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The action is expressed as 
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4 Numerical results 



The slave-boson equations (fT6j) and (fTUj) for a trapped condensate without and with vortex, 
respectively, can be solved numerically by means of a relaxation algorithm. The full 
solutions are compared to previous results, which were calculated from the Thomas-Fermi 
approximation (TFA) [TS]. In the TFA it is assumed that the kinetic energy can be 
neglected in comparison to the potential energy [26]. Thus the kinetic term containing the 
spatial derivatives is neglected. This leads to a transcendental equation which is easier to 
solve than the full differential equation. 

The relaxation algorithm which we employ now is similar to that applied for the GP 
equation in Ref. [28J, scheme B. First, we introduce an artificial dynamics by writing Eq. 
(|19p in a time-dependent form as 
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where \<f>\ 2 = \<f>\ 2 /rj, rj = J \cp\ 2 d s r. With this substitution the function 4> will always be 
normalized to 1 and bounded to values not too big or too small to be operated numerically 
by the kinetic term. This provides more stability to the code. 
Eq. (|22p is discretized in the split-step form 
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where the subscripts containing n correspond to the n th time step. Ocn corresponds to 
the Crank-Nicolson evolution algorithm alternated in r± and z directions [29J. In the 
rj_ direction, Neumann boundary conditions ( ^ T Q r ~°' z ^ = 0) were taken at the origin. 
The integrals of x an d Z can be evaluated efficiently at each time step by Gauss-Hermite 
quadrature with good convergence using only nine points. The algorithm starts with an 
initial ansatz distribution loaded to <p. After evolving for sufficient long time the time 
derivative Eq. (|22p vanishes, <j) n converges to the stationary solution and r] n converge to 
the norm of <p. 

In Fig|T]we show the condensate density profile of a trapped Bose gas at rest. As the 
chemical potential increases from (a) to (c), the total density of the Bose gas increases. It 
should be kept in mind that the total density is always monotonically decreasing from the 
center of the trap (cf. Ref. Due to increased interaction at higher total densities, 

depletation of the condensate reduces the condensate density at the center of the trap. As 
we can see in Fig{TJa), at small chemical potential (i.e. at low total density) the condensate 
density is monotonically decreasing from the center, but at higher values of the chemical 
potential depletion suppresses the condensate. This effect is stronger in regions with higher 
total densities (i.e. near the center of the trap), such that at sufficiently high total densities 
the condensate is entirely destroyed at the center (cf. Fig(Uc)). 

We also compare the full solution of Eq. (j!6f) and the TFA solution of Ref. [15] in 
Fig{TJ Although the difference between these solutions is small, the full solution shows that 
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the effect of the kinetic term is to smooth the abrupt changes of the TFA. Otherwise, the 
TFA is a very good approximation in the regime considered here. Moreover, the depletion 
of the condensate density at the center shows a perfect agreement between the full solution 
and the TFA in this region. 

In order to compare the slave-boson (SB) approach with resuls for the dilute Bose gas 
in the existing literature, we first evaluate the condensate density profile from the Gross- 
Pitaevskii (GP) equation (jlOp without (m = 0) and with (m = 1) vortex for the TFA and 
for the full kinetic GP equation. The results are shown in Figj2j Without a vortex we 
observe the smoothing of the kinetic term near the discontinuity around r = 5. In the 
vortex case, on the other hand, we find the formation of a core in the TFA. This core is 
also smoothed by the kinetic term in the full solution. 

Next, we evaluate the condensate density profile from the slave-boson approach, using 
Eq. ([?]), for corresponding parameter values. The results are presented in Figj3j The 
parameters where chosen so that the number of condensed atoms roughly agree in FigsJH 
and [3l 

These results show a general similarity of the GP and SB profiles. The cases of the 
GP and the SB calculations without vortex indicate that there is a noticeable depletion 
of the central condensate density in the SB approach. There is practically no difference 
between the TFA and the full solution. In the vortex case, the SB calculation shows that 
the core is also smooth when the kinetic term is added, as shown in the inset. Therefore 
the depletion in the center is overestimated by the TFA. This indicates that the formation 
of a vortex is less favorable than predicted by the TFA. Comparing the SB and the GP 
calculation in the vortex case, we also observe that the maximum of the SB profile is lower 
than that of the GP equation. Thus the depletion is more effective in the SB approach 
than in the GP equation. These results are in agreement with the general idea that the SB 
approach takes into account higher order interaction terms, leading to a more pronounced 
depletion at the center. 

The depletion at the center of the condensate should be observed experimentally by 
measuring the critical angular velocity of a stable vortex formation. In order to include 
strong interaction, the measurement should be performed near a Feshbach resonance, 
where the scattering length is very large. The critical angular velocity can be evaluated 
in the SB approach by following the procedure outlined previously in Ref. [9]. A vortex 
in a rotating condensate is stable if the action, given in Eq. ([3]), is smaller than the action 
of the condensate without vortex 

S m=1 (n) < S m=0 . (24) 

Solving numerically Eq. (j!9f) and using this stability criterion, we find the critical angular 
velocity f2 c for a vortex formation. In Fig. 2] we compile these results by plotting the 
critical angular velocity as a function of the number of condensed atoms. The TFA shows 
a lower critical angular velocity, both in the SB and in the GP approach, than the full 
numerical solutions. This is in qualitative agreement with previous results for the GP 
equation [30]. It can be simply understood in terms of our analysis above: Adding the 
kinetic term smooths the core generated by TFA. Thus the kinetic term adds condensed 
particles to the center. This effect suppresses the vortex formation. A comparison indicates 
that for a small number of atoms, i.e. low density, the SB prediction is practically the 
same as that of the GP equation. However, as we increase the number of atoms, the 
SB approach significantly deviates from the GP equation, showing lower critical angular 
velocity for a vortex formation. This can be explained with the profile analysis described 
above: the SB approach causes a stronger depletion than the GP equation. This difference 
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is even enhanced as the total density increases. As a consequence, a stronger depletion 
will favor a vortex formation and thus to a lower critical angular velocity. 

5 Conclusion 

We employed the relaxation method to solve the slave-boson model. With this method we 
obtained results beyond the Thomas-Fermi approximation by including the kinetic term 
of the field equation. We noticed that the depletion of the condensate is overestimated by 
the Thomas-Fermi approximation. Therefore, in this approximation the critical angular 
velocities are smaller than in the full solution. We also compared the full slave-boson model 
with the full Gross-Pitaevskii equation in the limit of low temperatures. In this limit, the 
slave-boson condensate profiles are very similar to the Gross-Pitaevskii condensate profiles 
at low densities. As the density is increased, there is a significant enhancement of the 
depletion in the slave-boson approach as compared to the Gross-Pitaevskii equation. This 
causes the critical angular velocity of vortex formation to be lower than expected from 
the Gross-Pitaevskii equation. The latter results corroborate previous estimates within 
the Thomas- Fermi approximation done in Ref.[l5]. All these effects should be observable 
in a Bose gas near a Feshbach resonance, where the scattering length is comparable with 
the mean distance of atoms. 
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Figure 1: Radial profile of the condensate density: Thomas-Fermi approximation results 
(dashed line) and the solution of the full differential equation of slave-boson approach 
(solid line) with /3' = 1 for the case without vortex (m = 0). Thomas- Fermi results were 
previously published in ref. [15]. (a) Chemical potential fjf = 0, (b) // = 1, (c) // = 2 
(units of J). The atomic density uq is given by equation (jUj) . r is the distance from the 
center of the trap and dh is the harmonic oscillator length. Labels are in dimensionless 
units, like in all the followig plots. (It should be noticed that the radial profile of the total 
density of particles is monotonous with its maximum at r = 




Figure 2: Condensate density profiles obtained from the Gross-Pitaesvkii (GP) equation 
using the Thomas-Fermi approximation (TFA) and full numerical solution for the cases 
without vortex (m = 0) and with vortex (m = 1). The parameters are the chemical 
potential \xqp = 0.5 and the total number of atoms N = 9095. r is the distance from the 
center of the trap and dho is the harmonic oscillator length. 
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Figure 3: Condensate density profiles obtained from slave-boson equation using Thomas- 
Fermi approximation (TFA) and full numerical solution without vortex (m = 0) and with 
vortex (m = 1). The parameters are (5' = 10, chemical potential [i = —4. 85 J, and number 
of condensed particles iVo = 9112. 
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Figure 4: Critical angular velocity calculated in Thomas-Fermi approximation and full cal- 
culations with kinetic term for the Gross-Pitaevskii (GP) and Slave-Boson (SB) equations. 
iVo is the number of condensed atoms. For the SB equation we used (3' = 10. 
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